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I. INTRODUCTION 



Structure functions in deep-inelastic scattering (DIS) and their scale evolution are closely 
related to the origins of quantum chromodynamics (QCD). DIS processes have played and 
still play a very important role for our understanding of QCD and nucleon structure In 
fact, DIS structure functions have been the subject of detailed theoretical and experimental 
investigations. Today, with high-precision data from the electron proton collider, HERA, 
and in view of the outstanding importance of hard scattering processes at proton(anti)proton 
colliders like the TEVATRON and the forthcoming Large Hadron Collider (LHC) at CERN, 
a quantitative understanding of deep-inelastic processes is indispensable. 

To predict the rates of the various processes, a set of universal parton distribution func- 
tions (PDF's) is required. On the other hand all calculations of high energy processes with 
initial hadrons, whether within the standard model or exploring new physics, require PDF's 
as an essential input. The reliability of these calculations, which underpins both future theo- 
retical and experimental progress, depends on understanding the uncertainties of the PDF's. 
These distribution functions can be determined by QCD global fits to all the available DIS 
and related hard-scattering data. The QCD fits can be performed at leading order (LO), 
next-to-leading order (NLO), next-to-next-to-leading order (N 2 LO) in the strong coupling 

The assessment of PDF's, their uncertainties and extrapolation to the kinematics relevant 
for future colliders such as the LHC have been an important challenge to high energy physics 
in recent years. Over the last couple of years there has been a considerable improvement 
in the precision, and in the kinematic range of the experimental measurements for many of 
these processes, as well as new types of data becoming available. In addition, there have 
been valuable theoretical developments, which increase the reliability of the global analysis. 
It is therefore timely, particularly in view of the forthcoming experiments at the LHC at 
CERN, to perform new global analysis which incorporate all of these improvements. A lot 
of efforts and challenges have been done to obtain PDF's for the LHC [2] which take into 
account the higher order corrections js-5]. 

For quantitatively reliable predictions of DIS and hard hadronic scattering processes, 
perturbative QCD corrections at the N 2 LO and the next-to-next-to-next-to-leading order 
(N 3 LO) need to be taken into account. Based on our experience obtained in a series of LO, 
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NLO and N 2 LO analysis [6j of the nonsinglet parton distribution functions, here we extend 
our work to N 3 LO accuracy in perturbative QCD. 

In this work this important problem is studied with the help of the method of the structure 
function reconstruction over their Mellin moments, which is based on the expansion of the 



structure function in terms o: 



for different QCD analyses 



case in Refs. 



7 



25| and [26|-|30| 



Jacobi polynomials. This method was developed and applied 



■24 ] . The same method has also been applied in the polarized 



In the present paper we perform a QCD analysis of the flavor nonsinglet unpolarized deep- 
inelastic charged e(/i)p and e{fi)d world data 3ll435l| at N 3 LO and derived parameterizations 
of valence quark distributions xu v (x, Q 2 ) and xd v (x, Q 2 ) at a starting scale Ql together with 
the QCD-scale Aq CD by using the Jacobi polynomial expansions. We have therefore used the 
3-loop splitting functions and Pade approximations |36l440| for the evolution of nonsinglet 
quark distributions of hadrons. 

Previous 3 loop QCD analysis were mainly performed as combined singlet and nonsin- 
glet analysis 4jJ, |42j, partly based on preliminary, approximative expression of the 3-loop 
splitting functions. Other analyses were carried out for fixed moments only in the singlet 



and nonsinglet case analyzing neutrino data 
were published in 



43 



Refs. Q, [|7, |48| . The results of 4-loop QCD analysis are also reported in |48|, |49j. The 



4451] . First results of the nonsinglet analysis 



46|| . Very recently a 3-loop nonsinglet analysis was also carried out in 



results of the present work are based on the Jacobi polynomials expansion of the nonsinglet 
structure function. 

The plan of the paper is to recall the theoretical formalism of the QCD analysis for 
calculating nonsinglet sector of proton structure function F 2 in Mellin-iV space in Sec. II. 
Section III explains the Pade approximations and 4-loop anomalous dimensions. A descrip- 
tion of the Jacobi polynomials and procedure of the QCD fit of F 2 data are illustrated in 
Sec. IV. The numerical results are illustrated in Sec. V before we summarize our findings in 
Sec. VI. 
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II. THEORETICAL FORMALISM OF THE QCD ANALYSIS 



In the common MS factorization scheme the relevant F 2 structure function as extracted 
from the DIS ep process can be written as 



50H53| 



x~ 1 F 2 (x,Q 2 



= x~ l (F 2)NS (x, Q 2 ) + F 2)S {x, Q 2 ) + F 2t9 (x, Q 2 )) 

= C 2 ,ns(x, Q 2 ) ® Qns(x, Q 2 ) 

+ < e 2 > C 2tS {x,Q 2 ) ®q s {x,Q 2 ) 

+ < e 2 > C 2>s (x, Q 2 ) ® g(x, Q 2 ) , 



(1) 



here and g represent the quarks and gluons distributions respectively, ^ns stands for the 
usual flavor nonsinglet combination and qs stand for the flavor- singlet quark distribution, 
qs = $2i=i(<ft + Qi)- Also, rij denotes the number of effectively massless flavors. < e 2 > 
represents the average squared charge, and £g> denotes the Mellin convolution which turns 
into a simple multiplication in iV-space. 

The perturbative expansion of the coefficient functions can be written as 



n=0 ^ ' 



X 



(2) 



In LO, C 2 ^ s (x) = 5(x), C^p S (x) = C 2 g(x) = C^p^x) = and the singlet-quark coefficient 



(i) 



(n) 
2 



function is decomposed into the nonsinglet and pure singlet contribution, C 
^ns + ^2PS- ^he coefficient functions C^f up to N 3 LO have been given in 54j 



a 



(n) 
2,S 



The nonsinglet structure function F 2jN s(x,Q 2 ) U P to N 3 LO and for three active (light) 
flavors has the representation 



x 1 F 2i ns(x,Q 2 



r (0) , r (l) , „2 r (2)+ 3^(3)+ 

°2,g 1^ u s L '2,NS 1^ "s°2,NS ' u, s°2,NS 

The flavor-singlet and gluon contributions in Eq. ([1]) reads 
x~ l F 2 , s {x,Q 2 ) = 





' 1 






.18 q 





C 2 J + a s C? I + a 2 C^ „ + Cr ^ 



, 's w 2,g 



's w 2,<j 



S(x,Q 2 



x _1 F2, s (x, Q 2 ) 
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+ n 2 r {2) + fl 3 r (3) ' 

a s L/ 2,g + a s L/ 2,g + CL s Ly 2,g 



The symbol <8> denotes the Mellin convolution 

[A<g>£](x) = 



1 ,1 

d^i / rfx 2 5(x — Xix 2 ) A(xi)B(x 2 ) . 
(I Jo 



(^,Q 2 ) • (3) 

(4) 
(5) 

(6) 



In Eq. q% = u+u — (d+d) = u v —d v and q£ = u+u+d+d—2(s + s) = u v +d v +2u+2d—As, 
where s — s. Also in Eq. (@|, E(x, Q 2 ) = S g=u ,d,s(? + q) — u v + d v + 2u + 2d + 2s. Notice 
that in the above equations a s = a s {Q^) = a s (Q 2 ) /An denotes the strong coupling constant 



and Cij are the Wilson coefficients 



54]. 



The combinations of parton densities in the nonsinglet regime and the valence region 



x > 0.3 for Fl in LO is 



x 



+ I + 

?NS,8 + g ^NS,3 



(x,Q 2 ) + -Z(x,Q 2 ) 



(7) 



where since sea quarks can be neglected 

in the region x > 0.3. So in the x-space we have 



^-xu v (x,Q 2 ) + \ xd v (x,Q 2 ) 

y y 



Fr?(x, Q 2 ) — f g a; 5ns,8 + g ^ ^ns.sJ ( x ' *3 ) 
In the above region the combinations of parton densities for F 2 are also given by 
F^(x,Q 2 ) = f^x?Ns, 8 ) ( X >Q 2 ) = Jg x ( u v + d v)(x,Q 2 ) , 



(8) 



(9) 



where g^s3 = u v — d v and F 2 d = (Ff + F£)/2 if we ignore the nuclear effects here. It is 
important to stress that the shadowing effect as a nuclear effect may affect our analysis. The 
shadowing effect 



55l . [56| arising from the gluon recombination and in the small- a; region 
the competitive mechanism of nuclear shadowing takes place. It also depends on the size of 
the nucleons. According to this effect we have F 2 = (F 2 + F 2 )/2 + 5F 2 . To obtain the 5F 2 
we need to know the generalized vector meson dominance (VMD) and parton mechanism at 
low and large values of Q 2 respectively. We found that the value of SF 2 is important but 
in low values of x. For example this correction value at Q 2 =10 GeV 2 and for x > 0.1 is too 
small (~ 10~ 4 ). So in the valence region of this analysis, this effect is negligible in large x 
and we can use the F 2 = (F 2 + F 2 )/2 approximately. 

In the region x < 0.3 for the difference of the proton and deuteron data we use 

F 2 NS (x,Q 2 ) = 2(F?-F 2 d )(x,Q 2 ) 

1 1 2 - 

= g ^s^Q 2 ) = 3 x ( u v ~ d v )(x,Q 2 ) + - x(u - d)(x,Q 2 ) , (10) 

d v + 2{u — d) since sea quarks cannot be neglected for x smaller 



where now q^s,3 
than about 0.3. 



u. 
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The first clear evidence for the flavor asymmetry combination of light parton distributions 



x(d — u) in nature came from the analysis of NMC at CERN to 



rule {57J . In our calculation we supposed the d — u distribution |47l . 



study o 



48 



the Gottfried sum 



58 



59) 



x(d — u)(x, Ql) = 1.195x 



1.24/ 



X 



,9.10, 



1 + 14.05x - 45.52x 2 



(11) 



at Ql = 4 GeV 2 which gives a good description of t 



601 ] . In this analysis, like other analyses 



ay 



47 



48 



re Drell-Yan dimuon production data 



58 



591 ) , we used the above distribution 



for considering the symmetry breaking of sea quarks. Although, in fact, this parametrization 
plays a marginal role in our analysis, in order to find the impact effect of this distribution, 
which is essentially used in the paper, it is desirable to study the QCD fits by varying 
this distribution with another asymmetry sea quark distribution which is derived in other 
analyses. In Sec. VI we will discuss our outputs when we change the above sea distribution. 

Now these results in the physical region < x < 1 can transform to Mellin-iV space by 
using the Mellin transform to obtain the moments of the structure function as ^F£, 



F*(N,Q 2 ) = M[F 2 fc , N] = I dx x^^-F^x, Q' 2 ) 

'0 x 



1 



(12) 



here k denotes the three above cases, i.e. k = p,d, NS. One of the advantages of Mellin- 
space calculations is the fact that the Mellin transform of a convolution of functions in 
Eqs. (!3|4f5l) reduces to a simple product 



M[A ®B,N] = M[A, N)M[B, N] = A(N)B(N) (13) 
By using the solution of the nonsinglet evolution equation for the parton densities to 4— 
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loop order, the nonsinglet structure functions are given by 48 1 

F 2 k (N, Q 2 ) = (l + a. Ci% s (N) + a 2 s Cgg(JV) + a 3 s Cg s (iV)) F 2 k (N, Q 2 ) 



x 



-Po(N)/(3 



1 - — (a s - a ) 

Po 



P+(n) - ^p (iv) 

Po 



1 



2po 
1 

Wo 
1 



a ) 



+ 



3po 
01 



P+(AT) - ^lp+(AT) 
Po 

P+(N) - ^P (iV) 
Po 



P 3 +(iV) - ^P 2 + (iV) + 
Po 



Pi 



ft 

PI 



Po 



Po 



Po{N) 



ao) (al-a 2 s ) (pt{N)-^h{N)j 



x 



h(N) - |A(*) - (I 



Po 



Po(iV) 



6/3 



«oJ 



01 



P+(iV) - ^P (W) 
Po 



(14) 



Here a s (= a s /47r) and a denotes the strong coupling constant in the scale of Q 2 and Ql 
respectively, k = p,d and NS also denotes the three above cases, i.e. proton, deuteron 
and nonsinglet structure function. C^ S (N) are the nonsinglet Wilson coefficients in 0{a™) 
which can be found in 
loop splitting functions. 



62 ] and P m denote also the Mellin transforms of the (m + 1)- 



III. PADE APPROXIMATIONS AND 4-LOOP ANOMALOUS DIMENSIONS 



In spite of the unknown 4-loop anomalous dimensions, one can obtain the nonsinglet 
parton distributions and Aq C d by estimating uncalculated fourth-order corrections to the 
nonsinglet anomalous dimension. On the other hand the 3-loop Wilson coefficients are 
known 54j and now it is possible to know, which effect has the 4-loop anomalous dimension 



if compared to the Wilson coefficient. In this case the 4-loop anomalous dimension may be 
obtained from Pade approximations. 

Pade approximations have proved to be useful in many physical applications. Pade 
approximations may be used either to predict the next term in some perturbative series, 



called a Pade approximation prediction, or to estimate the sum of the entire series, called 
Pade summation. 

For this purpose we use the Pade approximations of the perturbative series, discussed in 



detail for QCD, e.g., in Refs. 



36 



Pade approximations 



39 



40j are rational functions 



chosen to equal the perturbative series to the order calculated: 



[N/M] 



a + a\X + ... + a^x 1 ^ 
1 + hx + ... + b M x M 



to the series 



S = S + Six + ... + Sm +m x 



M+M 



(15) 



(16) 



where we set 



[N/M] = S + 0(x 



.Af+M+U 



and write an equation for the coefficients of each power of x. To continue, let's go to Mellin-iV 
space. 

A generic QCD anomalous dimension expansion in term of a s then may be written in the 
form 



7(A0 = $>* + V°(A0. 



(17) 



1=0 



In Mellin-iV space and by using this approach we can replace 7 (iV) by a rational function 
in a s [541 ] . 



1 + a s qi(N) + . . . + a-^qM^N) 



(18) 



Here M > 1 and N + M = n, where n stands for the maximal order in a s at which the 
expansion coefficients 7 <n )(iV) have been determined from an exact calculation. The func- 
tions Pi(N) and qj{N) are determined from these known coefficients by expanding Eq. ([18]) 
in powers of a s . This expansion then also provides the [N/M] Pade approximate for the 
(n+l)-th order quantities ^ n+1 \ 

In this way it is easy to obtain the following results for M = N = 1 and for M = 0, N = 2 

7 (2)2 (iV) 



7 [1/11 (iV) = [1/1] (N) 
~[o/ 2](iV) _ [o/2] (AT) 



7 ( 1 )(A^) ' 

2 7 ( 1 )(AT) 7 ( 2 )(iV) ^(N) 



7 (°)(iV) 



7 (°) 2 (AT) 



(19) 



The strong coupling constant a s plays a more central role in the present paper to the 
evolution of parton densities. At N m LO the scale dependence of a s is given by 

da s 



dlnQ 2 



f3N m Lo(as) 



k=0 



,k+2 



Ik 



(20) 



The expansion coefficients /3 k of the /3-function of QCD are known up to k = 3, i.e., N 3 LO 
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64] 



A, = 11 - 2/3 n f , 

0! = 102 - 38/3 n f , 

(3 2 = 2857/2 - 5033/18 n f + 325/54 n) , 

/3 3 = 29243.0 - 6946.30 n f + 405.089 n) + 1093/729 n) , (21) 
here n/ stands for the number of effectively massless quark flavors and /3k denote the coeffi- 



cients of the usual four- dimensional MS beta function of QCD. In complete 4-loop approx 



imation and using the A-parametrization, the running coupling is given by 
1 1 



65 



66J: 



a s (Q 2 



+ 
+ 



P L A 
1 



bi In L[ 



[bl (ln 2 L A -lnL A -l) + b 2 ] 



3btb 2 In L A + — 



(22) 



5 1 
In 3 L A + -ln 2 L A + 21nL A -- 

where L A = ln(Q 2 /A 2 ), bk = fik/Po, and A is the QCD scale parameter. The first line of 
Eq. (1221) includes the 1- and the 2-loop coefficients, the second line is the 3-loop and the 
third line denotes the 4-loop correction. Equation ( 1221) solves the evolution equation ( 1201) 
only up to higher orders in 1/L A . The functional form of a s (Q 2 ), in 4-loop approximation 
and for 6 different values of A, is displayed in Fig. [TJ The slope and dependence on the 
actual value of A is especially pronounced at small Q 2 , while at large Q 2 both the energy 
dependence and the dependence on A becomes increasingly feeble. To be able to compare 
with other measurements of A we adopt the matching of flavor thresholds at Q 2 = m 2 and 
Q 2 = ml with m c = 1.5 GeV and = 4.5 GeV as described in [67|, l68 |. 



IV. JACOBI POLYNOMIALS AND THE PROCEDURE OF QCD FITS 



One of the simplest and fastest possibilities in the structure function reconstruction from 
the QCD predictions for its Mellin moments is Jacobi polynomials expansion. The Jacobi 

9 



polynomials are especially suitable for this purpose since they allow one to factor out an 
essential part of the x-dependence of structure function into the weight function 8j. 

According to this method, one can relate the F 2 structure function with its Mellin mo- 
ments 

■^V max 

F 2 "> N <™ (*,Q 2 ) = x^l-xT £ e^(x)^cf ) (a,/3)F 2 fc ( J + 2,Q 2 ), (23) 

71=0 j = 

where N max is the number of polynomials, k denotes the three cases, i.e. k — p,d, NS. 



Jacobi polynomials of order n 69j, Q" ,l3 (x), satisfy the orthogonality condition with the 
weight function w af5 = ar (1 — x) a 

I dx w a ^e^(x)ef' P (x) = 5 k>l . (24) 
Jo 

In the above, c^\a,f3) are the coefficients expressed through T-functions and satisfying the 
orthogonality relation in Eq. (j24j) and F 2 (j + 2,Q 2 ) are the moments determined in the 
previous section. N max , a and (3 have to be chosen so as to achieve the fastest convergence 
of the series on the right-hand side of Eq. (1231) and to reconstruct F 2 with the required 
accuracy. In our analysis we use N max = 9, a = 3.0 and (3 = 0.5. The same method has 



been applied to calculate the nonsinglet structure function xF 3 from their moments 



13- 



16| and for polarized structure function xg± [25l-l27l|. Obviously the Q 2 -dependence of the 



polarized structure function is defined by the Q 2 -dependence of the moments. 

The evolution equations allow one to calculate the Q 2 -dependence of the parton distri- 
butions provided at a certain reference point Q^. These distributions are usually parame- 
terized on the basis of plausible theoretical assumptions concerning their behavior near the 
end points x = 0,1. 

In the present analysis we choose the following parametrization for the valence quark 
densities in the input scale of Qq = 4 GeV 2 

xq v (x,Q 2 ) =Af q x a «(l-x) b o(l + c q y/x~ + d q x) , (25) 

where q = u,d and the normalization factors Af u and Afd are fixed by J Q u v dx = 2 and 
Jq 1 d v dx = 1, respectively. By QCD fits of the world data for F$' d , we can extract valence 
quark densities using the Jacobi polynomials method. For the nonsinglet QCD analysis 
presented in this paper we use the structure function data measured in charged lepton-proton 



10 



and deuteron deep-inelastic scattering. The experiments contributing to the statistics are 



BCDMS 



3l|, SLAC [32J, NMC 



33], HI [34| . and ZEUS |35|. In our QCD analysis we use 
three data samples : F%(x, Q 2 ), F${x, Q 2 ) in the nonsinglet regime and the valence quark 
region x > 0.3 and F^ s = 2(F 2 — -F 2 d ) in the region x < 0.3. 

The valence quark region may be parameterized by the nonsinglet combinations of parton 
distributions, which are expressed through the parton distributions of valence quarks. Only 
data with Q 2 > 4 GeV 2 were included in the analysis and a cut in the hadronic mass of 
W 2 = (- — 1) Q 2 + m^j > 12.5 GeV 2 was applied in order to widely eliminate higher twist 
(HT) effects from the data samples. After these cuts we are left with 762 data points, 322 
for F 2 P , 232 for F 2 d , and 208 for F* s . By considering the additional cuts on the BCDMS 
(y > 0.35) and on the NMC data(Q 2 > 8 GeV 2 ) the total number of data points available 
for the analysis reduce from 762 to 551, because we have 227 data points for F 2 , 159 for F%, 
and 165 for F^ s . 

For data used in the global analysis, most experiments combine various systematic errors 
into one effective error for each data point, along with the statistical error. In addition, the 
fully correlated normalization error of the experiment is usually specified separately. For 
this reason, it is natural to adopt the following definition for the effective \ 2 



7Q| 



Xgiobai = w nXn > { n labels the different experiments) 

n 

i - M n \ 2 ^ ( KF^ ta - F' h r r \ 2 



2 / x • ,v " 



For the n th experiment, F^f ", AF 2 rf f a , and F^ OT denote the data value, measurement 
uncertainty (statistical and systematic combined) and theoretical value for the i th data point. 
AJ\f n is the experimental normalization uncertainty and M n is an overall normalization factor 
for the data of experiment n. The factor w n is a possible weighting factor (with default value 
1). However, we allowed for a relative normalization shift M n between the different data 
sets within the normalization uncertainties AJ\f n quoted by the experiments. For example 
the normalization uncertainty of the NMC(combined) data is estimated to be 2.5%. The 
normalization shifts M n were fitted once and then kept fixed. 

Now the sums in Xgiobai run over a ^ data sets and in each data set over all data points. The 
minimization of the above x 2 value to determine the best parametrization of the unpolarized 



parton distributions is done using the program MINUIT 

11 



7l|. 



The one a error for the parton density xq v as given by Gaussian error propagation is 48] 



a{xq v {x)f = jrj^ (^) C0V (PiiPi) > ( 27 ) 

where the sum runs over all fitted parameters. The functions dxq v /dpi are the derivatives 
of xq v with respect to the fit parameter p i} and cov(pi,pj) are the elements of the covariance 
matrix. The derivatives dxq v /dpi can be calculated analytically at the input scale Qq. Their 
values at Q 2 are given by evolution which is performed in Mellin-iV space. 

Now we need to discuss the derivatives in Mellin-iV space a bit further. The Mellin-iV 
moment for complex values of N calculated at the input scale Ql for the parton density 
parameterized as in Eq. fl25|) is given by 



q v (N,a q ,bq,c g ,d q ) = J\f q M(n,a q ,b q ,c q ,d q ) , (28) 
with the normalization constant 

a 

j\f = z±i (29) 

q M(l, a q ,b q ,c q ,d q ) ' 

Here C Qv is the respective number of valence quarks, i.e. C Uv =2 and C^=l. In the above 
M(n, a q , b q , c q , d q ) is given by 

M(n, a qi b qi c qi d q ) = B[a q + n — 1, b q + 1] + c q B[a + n + 1/2, b + 1] + d u B[a q + n, b q + 1] , 

(30) 

where B[a,b] denotes the Euler beta function for complex arguments. The general form of 
the derivative of the Mellin moment q v with respect to the parameter p is given by 

^ = „,»,a +Af ,?M, (31) 
op op op 

In this analysis only the parameters a q and b q have been fitted for both the xu v and xd v 
parametrization while the other parameters involved are kept fixed after a first minimization 
in the MINUIT program, since their errors turned out to be rather large compared to the 
central values. Here we want to show the derivatives u v and d v parton densities with respect 
to parameter a q and b q . For example: 

<9M(n, a ) 

f(n,a q ) = ' q = B[a q +n - l,b q + l]{if)[a q + n - 1] - if)[a q + b q + n]) + 

c q B[a q +n- 1/2, b q + l](tp[a q + n- 1/2] - ip[a + b + n + 1/2]) + 

d q B[a q + n,b q + l](ip[a q + n] — ijj[a q + b q + n + 1]) , (32) 

12 



= B[a q +n-l,b q + l){ip[b q + 1] - ^P[a q + b q + n)) + 

c q B[a q + n- 1/2, b q + l}(ip[l + b q ] - tp[a q + b q + n + 1/2]) + 
d q B[a q + n,b q + l](ip[b q + 1] - 4>[a q + b q + n + 1]) , 




(33) 



and now we can reach the below derivatives for u v (N) and d v (N) with respect to parameters 
a q and b q 



also ip[n] = d In T{n)/dn is Euler's ^-function. 

To obtain the error calculation of the structure functions Ff , F$ , and F? s the relevant 
gradients of the PDF's in Mellin space have to be multiplied with the corresponding Wilson 
coefficients. This yields the errors as far as the QCD parameter A is fixed and regarded 
uncorrelated. The error calculation for a variable A is done numerically due to the nonlinear 
relation and required iterative treatment in the calculation of a s (Q 2 , A) {(], Q. 




V. RESULTS 

In the QCD analysis of the present paper we used three data sets: the structure functions 
F^ix, Q 2 ) and F^ix, Q 2 ) in the region of x > 0.3 and the combination of these structure 
functions F| fS (a;,Q 2 ) in the region of x < 0.3 . Notice that we take into account the cuts 
Q 2 > 4 GeV 2 , W 2 > 12.5 GeV 2 for our QCD fits to determine some unknown parameters. 
In Fig.(|2]) the proton, deuteron and nonsinglet data for Q 2 ), F^x, Q 2 ) and Ff s (x, Q 2 ) 
are shown in the nonsinglet regime and the valence quark region x > 0.3 indicating the above 
cuts by a vertical dashed line. The solid lines correspond to the N 3 LO QCD fit. Now, it is 
possible to take into account the target mass effects in our calculations. The perturbative 
form of the moments is derived under the assumption that the mass of the target hadron 
is zero (in the limit Q 2 — > oo). At intermediate and low Q 2 this assumption will begin to 
break down and the moments will be subject to potentially significant power corrections, of 
order O (mfj/Q 2 ), where is the mass of the nucleon. These are known as target mass 
corrections (TMCs) and when included, the moments of flavor nonsinglet structure function 



dg v (N, p) 
dp 



K(f(n,p)-f(l,p)M(n,p)/M(l,p)) , 



(34) 
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have the form 



43, LSI 

k 2 f 1 n-1 1 k 2 

JO x 



'0 

= ff ( B ,g») + 2^1(^) f J( B + 2,«") 

where higher powers than (m^/Q 2 ) 2 are negligible for the relevant x < 0.8 region. By 
inserting Eq. (I3"5"j) in Eq. (1231) we have 

ax 71 

F 2 k > N ™(x,Q 2 ) = x^l-x)* Q n^ x ) >< Ecf ) (a,/3)F 2 fc TMC ( J + 2,Q 2 ) , (36) 

n=0 j=0 

where -FItmcC/ + 2, Q 2 ) are the moments determined by Eq. fl35l) . In Fig.(j5J) the dashed 
lines correspond to the N 3 LO QCD fit adding target mass corrections. 

Despite the kinematic cuts (Q 2 > 4 GeV 2 , W 2 = (i - 1)Q 2 + m& > 12.5 GeV 2 ) 
used for our analysis, we also take into account higher twist corrections to -Ff (x, Q 2 ) and 
F^x^Q 2 ) in the kinematic region Q 2 > 4 GeV 2 , 4 < W 2 < 12.5 GeV 2 in order to learn 
whether nonperturbative effects may still contaminate our perturbative analysis. For this 
purpose we extrapolate the QCD fit results obtained for W 2 > 12.5 GeV 2 to the region 
Q 2 > 4 GeV 2 , 4 < W 2 < 12.5 GeV 2 and from the difference between data and theory, 
applying target mass corrections in addition. Now by considering higher twist correction 

FTix, Q 2 ) = TMC [if T (*, Q 2 )] ■ (l + , (37) 

the higher twist coefficient can be extract. Here the operation Otmc[--] denotes taking the 
target mass corrections of the twist-2 contributions to the respective structure function. 
The coefficients h(x, Q 2 ) are determined in bins of x and Q 2 and are then averaged over 
Q 2 . We extrapolate our QCD fits to the region 12.5 GeV 2 > W 2 > 4 GeV 2 in Fig.©. The 
dash-dotted lines in this figure correspond to the N 3 LO QCD fit adding target mass and 
higher twist corrections. There, at higher values of x a clear gap between the data and the 
QCD fit is seen. 



In Table $ we summarize the NLO, N 2 LO, and N 3 LO with using Pade [1/1] and [0/2] 
fit results for the parameters of the parton densities xu v (x, Qq), xd v (x, Qq) and Aqq^. The 
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NLO 


N 2 LO 


N 3 LO Pade [1/1] 


N 3 LO Pade [0/2] 


u v 


a u 
b u 

C-u 

d u 


0.7434 ± 0.009 
3.8907 ± 0.040 

0.1620 

1.2100 


0.7772 ± 0.009 
4.0034 ± 0.033 

0.1000 

1.1400 


0.79167 ± 0.0106 
4.02637 ± 0.0402 

0.0940 

1.1100 


0.79176 ± 0.0099 
4.02685 ± 0.0327 

0.0940 

1.1100 




ad 
b d 

Cd 

dd 


0.7369 ± 0.040 
3.5051 ± 0.225 

0.3899 

-1.3700 


0.7858 ± 0.043 
3.6336 ± 0.244 

0.1838 

-1.2152 


0.80927 ± 0.0621 
3.76847 ± 0.3499 

0.1399 

-1.1200 


0.80927 ± 0.0407 
3.76858 ± 0.2278 

0.1399 

-1.1200 


AJcd 4 , MeV 


263.8 ± 30 


239.9 ± 27 


241.44 ± 29 


241.45 ± 27 


X 2 /ndf 


523/546 = 0.9578 


506/546 = 0.9267 


491.07/546 = 0.8994 


491.12/546 = 0.8995 



TABLE I: Parameter values of the NLO, N 2 LO from Ref. [6] and N 3 LO nonsinglet QCD fit at Ql = 
4 GeV 2 for Pade [1/1] and Pade [0/2]. 

values without error have been fixed after a first minimization since the data do not constrain 
these parameters well enough. In this table we also compare the N 3 LO results with the NLO 
and N 2 LO results from Ref.[6|. The results show a good compatibility between Pade [1/1] 
and [0/2] approximations in 4-loop order. The resulted value of x^/ndf is 0.9578 at NLO, 
0.9267 at N 2 LO, and 0.8994 and 0.8995 for Pade [1/1] and [0/2] respectively at N 3 LO. Our 
results for the covariance matrix of the N 3 LO nonsinglet QCD fit for Pade [1/1] and [0/2] 
are presented in Table (ILT1). 

Figure ([H illustrates our fit results for x u v (x, Ql), xd v (x, Ql) at Ql = 4 GeV 2 up to N 3 LO 



and for Pade [1/1] wit 
results obtained from 



i correlated errors. In this figure our results for N 3 LO compared with 
6[ at LO, NLO, and N 2 LO QCD analysis. The shaded areas represent 
the fully correlated one a statistical error bands. 

In Fig. 01]) we show the evolution of the valence quark distributions xu v (x,Q 2 ) and 
xd v (x, Q 2 ) from Q 2 = 1 GeV 2 to Q 2 = 10 4 GeV 2 in the region x G [10" 4 , 1] at N 3 LO. In this 



48[ . With rising 



figure we also compared our results with the nonsinglet QCD analysis from 
values of Q 2 the distributions flatten at large values of x and rise at low values. 

Another way to test the N 3 LO fit results is comparison of low order moments of the 
distributions u v (x, Q 2 ), d v (x, Q 2 ), and u v (x,Q 2 ) — d v (x,Q 2 ). In Table [LTTl we present the 
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TV t ~w~ y — ^ t~\ i / r -1 / -t ~\ 

N 3 LO Pade[l/1] 


a u 


b u 




b d 


a N f =4 
QCD 




1.13xl0~ 4 










K 


O or , , 1 r\ — 4 

2.35x11) 


1.62x10 








ad 


i.uy x iu 


-i.oy x iu 


O.OO x ±u 






t>d 


1.0 / X 1U 


-o.o4fc X IU 


o ii vi n — 2 

Z.11X1U 






A (4) 
QCD 


1 71 v i n — 4 
1 . / 1X1U 


Q /in vi n— 4 

-o.4y x iu 


cn/ivi n — 4 

D.U4: X IU 


Z.Dl X IU 


o cKvi n~4 

O.OO X _LU 


N 3 LO Pade[0/2] 








b d 


a N f =4 
QCD 


a u 


0.98x10 4 












1 SQvin - 4 
l.OO X 1U 


±.U I x ±u 








ad 


-5.07xl(T 5 


-6.01 xl0~ 4 


1.66x10 3 






b d 


-l.llxlO" 4 


-3.30xl0~ 3 


8.58xl(T 3 


5.19xl0~ 2 




A (4) 


1.59xl0" 4 


-1.99xl0~ 4 


1.94xl0" 4 


8.07xl0" 4 


7.53xl0~ 4 



TABLE II: Our results for the covariance matrix of the N 3 L0 nonsinglet QCD fit for Pade [1/1] and 



[0/2] at Ql = 4 GeV 2 by using MINUIT 



lowest nontrivial moments of these distributions at Q 2 = Ql in N 3 LO and compare to the 
respective moments obtained for the parameterizations 48]. 

We should note that the unknown parameters are correlated and almost depend on the 
method of the QCD fits. We believe that the source of the small difference between the 
results of our analysis and reported results in [48j is the kind of the different method of the 
QCD analysis. We used the Jacobi polynomial method as an expansion method to do QCD 
fits but they used the exact inverse Mellin technique to obtain some unknown parameters. 
We also found that the results of Pade [1/1] and [0/2] in 4- loop level are almost the same. 

To perform higher twist QCD analysis of the nonsinglet world data in N 3 LO, we consider 
the Q 2 > 4 GeV 2 , 4 < W 2 < 12.5 GeV 2 cuts. The number of data points in the above range 
for proton and deuteron is 279 and 278, respectively. The extracted distributions for h(x) 
in N 3 LO are depicted in Fig. ([5]) for the nonsinglet case considering scattering off the proton 
and deuteron target. According to our results the coefficient h(x) grows towards large x.To 
compare, we also present the reported results of the early N 2 LO analysis [6] in Fig.(j5]). Also 
in this figure HT contributions have the tendency to decrease form N 2 LO to N 3 LO. This 

n 

effect was observed for the first time in the case of fits of F3 DIS vN data in uM and then 
studied in more detail in 15 
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/ 


N 


BBG [48] 


N 3 LO Pade[l/1] 


N 3 LO Pade[0/2] 


u v 


2 


0.3006 ± 0.0031 


0.30757 ± 0.0026 


0.30806 ± 0.0028 




3 


0.0877 ± 0.0012 


0.08771 ± 0.0011 


0.08781 ± 0.0012 




4 


0.0335 ± 0.0006 


0.03320 ± 0.0006 


0.03323 ± 0.0006 


d v 


2 


0.1252 ± 0.0027 


0.12450 ± 0.0024 


0.12495 ± 0.0025 




3 


0.0318 ± 0.0009 


0.03040 ± 0.0008 


0.03012 ± 0.0008 




A 

4 


U.UIUO ± U.UUU4 


u.uuyyz ± U.UUU4 


u.uuyyo ± u.uuuo 




2 


0.1754 ±0.0041 


0.18305 ± 0.0036 


0.18310 ± 0.0038 




3 


0.0559 ± 0.0015 


0.05767 ± 0.0013 


0.05769 ± 0.0014 




4 


0.0229 ± 0.0007 


0.02329 db 0.0007 


0.02329 ± 0.0007 



TABLE III: Comparison of low order moments from our nonsinglet N 3 LO QCD analysis at Qq = 4 GeV 2 



with the N 3 LO analysis from Ref. 



48] 





a 


b 


c 


Proton 


1.015 


3.928 


- 0.193 


Deuteron 


4.481 


7.759 


- 0.064 



TABLE IV: Our results for h(x) function according to Eq. ([38]) and for N 3 LO. 



This similar effect was also observed in the fits of F 2 charge lepton-nucleon DIS data 



[471, 



4s[ HQ- 

In Ref. j^l, the functional form for h(x) is chosen by 



h(x) = a 



1 — x 



(38) 



and it is possible to compare h(x) results in N 2 LO and N 3 LO . In Table HVl we present our 
results for a, b, c in the above equation. 

As seen from Fig. (jSJ) h(x) is widely independent of the target comparing the results for 
deeply inelastic scattering off protons and deuterons. 



VI. DISCUSSION 

A study [6] of the available world data on deep-inelastic lepton-proton and lepton- 
deuteron scattering provided a determination of the valence quark parton densities and 
a s in wide ranges of the Bjorken scaling variable x and Q 2 up to 3-loop. In the nonsinglet 
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case, where heavy flavor effects are negligibly small, the analysis can be extended to 4-loop 
level, i.e. to QCD in N 3 LO perturbative expansion. 

The analysis was performed using the Jacobi polynomials method to determine the pa- 
rameters of the problem in a fit to the data. A new aspect in comparison with previous 
analysis is that we determine the parton densities and the QCD scale up to N 3 LO by us- 
ing the Jacobi polynomial expansion method and using Pade approximations. The benefit 
of this approach is the possibility to determine nonsinglet parton distributions analytically 
and not numerically. In Ref. [75| we arrange the MATHEMATICA program to extract 
xu v (x, Q 2 ) and xd v (x,Q 2 ) up to the 4-loops. 



In this analysis we adopt the d — u distribution at Qq 



> 2 - 4 GeV 2 from Refs. 
which gives a good description of the Drell-Yan dimuon production data 



47 



48 



58 



59j, 



60] . The nonsinglet 



regime is manifesting itself at x > 0.1 as the rule. In this regime, when we changed the sea 
distribution from the other groups, the value of x 2 , valence distributions, A and a s varied, 
but only slightly. For example, we used the d — u distribution from {3], 76-78] and we found 
that the value of x 2 varies by about 3% and A by about l%-2%. 

In the QCD analysis we parameterized the strong coupling constant a s in terms of four 
massless flavors determining Aqcd- Up to N 3 LO results fitting the data, are 



A (4) 
QCD 

A (4) 
^QCD 

A (4) 
QCD 



A 



(4) 
QCD 



213.2 ± 28 MeV, LO, 

263.8 ± 30 MeV, NLO, 

239.9 ±27 MeV, N 2 L0, 
241.4 ±29 MeV, N 3 L0. 



(39) 



These results can be expressed in terms of a s (M|): 



«s(Mf) 
«s(Mf) 
« s (Mf) 
a s (M 2 z ) 



0.1281 ±0.0028, L0, 

0.1149 ±0.0021, NLO, 

0.1131 ±0.0019, N 2 L0, 

0.1139 ±0.0020, N 3 L0. 



(40) 



Note that in above results we use the matching between rif and flavor couplings calcu- 
lated in Ref. [661]. We adopt this prescription to be able to compare our results with other 
measurement of Aq CD . 
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The a s (M§) values can be compared with results from other QCD analysis of inclusive 
deep-inelastic scattering data in N 2 LO 



A02 [42]: a s (Mf)=0.1143 ±0.0014 

GRS [47]: a s (Mf)=0.111 

MRST03 [41]: a s (Mf) = 0.1153 ±0.0020 

SYOl(ep) [43]: a a (Mf )=0.1166 ±0.0013 

SY01(z/N) [43]: a s (M|)=0.1153 ±0.0063 

A06 [79]: a s (Mf )=0.1128 ± 0.0015 

BBG [48]: a s (Mf)=0.1134 + °'°° 19 

-0.0021 

BM07 [80]: a s (Mf )=0.1189 ±0.0019 
KPS00(z/Y) [15]: a s (Mf )=0.118 ± 0.002 (stat)± 0.005 (syst) 

± 0.003 (theory) 
KPS03(WV) [16]: a s (Mf)=0.119 ± 0.002 (stat)± 0.005 (syst) 

± 0.002 (threshold) 1°;°°* (scale) 
KT08 [6J: a s (Mf)=0.1131 ±0.0019 

The N 3 LO values of a s (M|) can also be compared with results from other QCD analysis 



BBG [48]: a s (Mf)=0.1134 + °'°° 19 

-0.0021 



and with the value of the world average 0.1189 ± 0.0010 8l| and also the current world 
average 

a s (M 2 z ) = 0.1184 ±0.0007 , (41) 

which has been extracted in 



82] very recently. It seems that our results confirm that the 
value of a s (M§) from DIS turns out to be sizably below the world average. In this case, 
it would be useful to find out which set of data is mainly responsible for the low value of 
a s (M§). We will try to see which subset makes a s (M|) particularly small in a future work. 

We hope our results of QCD analysis of structure functions in terms of Jacobi polynomials 
could be able to describe more complicated hadron structure functions. We also hope to be 
able to consider massive quark contributions by using the structure function expansion in 
terms of the Jacobi polynomials. 
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FIG. 1: The strong running of a s (Q 2 ), according to Eq.[22l in 4- loop approximation and for different 
values of A. 
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FIG. 2: The structure functions F% , F$, and FJ? s as a function of Q 2 in intervals of x. Shown are the 
Pade [1/1] QCD fits in N 3 LO (solid line) and the contributions from target mass corrections (dashed 
line) and higher twist (dash-dotted line). The vertical dashed line indicates the regions with W 2 > 12.5 
GeV 2 . 
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FIG. 3: The parton densities xu v and o;d„ up to 4-loop (Pade [1/1]) at the input scale Qq = 4.0 GeV 2 
(solid line) compared with results obtained from N 2 LO analysis (dashed- line), NLO analysis (dash- 
dotted line), and LO analysis(dash-dott-dotted line) [6]. The shaded areas represent the fully correlated 
one a statistical error bands. 
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FIG. 4: The parton densities xu v and xd v at N 3 L0 evolved up to Q 2 = 10000 GeV 2 (solid lines) 



compared with results obtained by BBG (dashed line) 
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FIG. 5: The higher twist coefficient h(x) for the proton and deuteron data as a function of x in 
N 3 LO (solid line) compared with results obtained by N 2 LO (dashed line) [f| . 
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